Evaluation of genomic offset predictions in common gardens

Author

Juliette Archambeau

Published

January 18, 2024

1 Introduction

In this report, we evaluate whether genomic offset (GO) predictions from the different GEA methods (LFMM, GF, GDM and RDA) and SNP sets (control, candidate and all SNPs for LFMM) are good predictors of fitness proxies (height and mortality data) from five clonal common gardens located in Spain (Asturias, Madrid, Cáceres), Portugal (Fundão) and France (Pierroton). For that, we predicted the genomic offset of each population when transplanted in the environment of the common garden (so instead of predicting the genomic offset into future climates, we predict it in the new environment of the common garden, i.e. space-for-time approach). The rational is that populations with higher genomic offset in a given common garden are expected to have lower relative fitness in the common garden.

We also calculate the climatic transfer distances (CTD) between the climate-of-origin of the populations and the climate of the common garden (i.e. the absolute value of the difference between the climate of origin of the populations and the climate in the common garden between the tree planting date and the measurement date). We then estimate the association between tree height and mortality and the climatic transfer distances. The climatic transfer distances are calculated for the same set of climatic variables as the one used to calculate the genomic offset (one CTD per climatic variable).

2 Data

2.1 Genomic offset predictions

We load the genomic offset predictions estimated from the different GEAs.

Code
df <- lapply(c("control","cand"), function(x){

list_snps <- list()

list_snps[[x]]$GDM <- readRDS(file=here("outputs/GDM/go_predictions.rds"))[[x]][["go_cg"]]
list_snps[[x]]$GF <- readRDS(file=here("outputs/GF/go_predictions.rds"))[[x]][["go_cg"]]
list_snps[[x]]$RDA <- readRDS(file=here("outputs/RDA/go_predictions.rds"))[[x]][["go_cg"]]
list_snps[[x]]$LFMM <- readRDS(file=here("outputs/LFMM/go_predictions_snpsets.rds"))[[x]][["go_cg"]]


list_snps <- list_snps[[x]] %>% 
  bind_rows(.id="method_type") %>% 
  mutate(method_input = x)

return(list_snps)
}) %>% bind_rows()


df <- readRDS(file=here("outputs/LFMM/go_predictions_allsnps.rds"))[["go_cg"]] %>% 
   mutate(method_type = "LFMM",
          method_input = "all") %>% 
  bind_rows(df) %>% 
  pivot_longer(cols=c(readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds"))[["cg"]]),names_to="cg",values_to="varX") %>% 
  mutate(method = paste0(method_type, "_",method_input))

2.2 Climatic transfer distances

We calculate the climatic transfer distances.

Code
# selected climatic variables
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))

# climatic data in the common gardens (btw planting and measurement dates)
cg_clim <-readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% 
  dplyr::select(cg,any_of(clim_var)) 

# Loading point estimate climatic data
adj <- "noADJ"  # not adjusted for elevation
ref_period <- "ref_1901_1950" # reference period 1901-1950
clim_past <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,any_of(clim_var))


df <- lapply(cg_clim[["cg"]], function(x){
  
for(var in clim_var){
  
  clim_past[[var]] <- ( clim_past[[var]] - cg_clim %>% filter(cg == x) %>%  pull(var) ) %>% abs()
} 
  
return(clim_past)
}) %>%  
  setNames(cg_clim[["cg"]]) %>% 
  bind_rows(.id="cg") %>% 
  pivot_longer(cols=any_of(clim_var),names_to="method_input",values_to="varX") %>% 
  mutate(method_type="CTD",
         method=paste0("CTD_",method_input)) %>% 
  bind_rows(df)

2.3 Phenotypic data

We load the survival and mortality data from the five common gardens.

Code
pheno_data <- readRDS(file=here("data/CommonGardenData/PhenoDataNovember2019_AnnualTraits_UpdatedSept2021_AllSites.rds")) %>% dplyr::rename(pop=prov)

no_nas <- sapply(pheno_data, function(x) length(x)-sum(is.na(x)))

list_pheno <- list()

list_pheno$`Asturias, Spain (37 months)` <- table(pheno_data$site,pheno_data$AST_survmar14)["asturias",]
list_pheno$`Bordeaux, France (85 months)` <- table(pheno_data$site,pheno_data$BDX_surv18)["bordeaux",]
list_pheno$`Cáceres, Spain (8 months)` <- table(pheno_data$site,pheno_data$CAC_survdec11)["caceres",]
list_pheno$`Madrid, Spain (13 months)` <- table(pheno_data$site,pheno_data$MAD_survdec11)["madrid",]
list_pheno$`Fundão, Portugal (27 months)` <- table(pheno_data$site,pheno_data$POR_survmay13)["portugal",]

df_exp_design <- list_pheno %>% 
  bind_rows(.id="cg") %>% 
  setNames(c("Common garden (tree age)",
             "Number of dead trees",
             "Number of trees alive")) %>% 
  mutate("Number of height measurements"=c(no_nas[["AST_htmar14"]],
                                           no_nas[["BDX_htnov18"]],
                                           no_nas[["CAC_htdec11"]],
                                           no_nas[["MAD_htdec11"]],
                                           no_nas[["POR_htmay13"]]))

saveRDS(df_exp_design, file = here("tables/ExpDesign_CG.rds"))

df_exp_design  %>% kable_mydf
Common garden (tree age) Number of dead trees Number of trees alive Number of height measurements
Asturias, Spain (37 months) 206 3976 3976
Bordeaux, France (85 months) 119 3225 3217
Cáceres, Spain (8 months) 3818 366 340
Madrid, Spain (13 months) 3138 1046 1046
Fundão, Portugal (27 months) 1517 2666 2665

Height and survival measurements in each common garden:

  • Asturias (Spain) in March 2014 when the trees were 37 month-old (trees were planted in February 2011).

  • Bordeaux (France) in November 2018 when the trees were 85 month-old (trees were planted in October 2011).

  • Cáceres (Spain) in December 2011 when the trees were 8 month-old (trees were planted in April 2011). Note that for this common garden, we calculate the bioclimatic variables for the entire year 2011 (instead of calculating the variables only for months between April and December). Indeed, the calculation of the annual bioclimatic variables will be wrong if we do not account for some months, e.g. the mean annual temperature will be higher than expected because we do not account for some winter months.

  • Madrid (Spain) in December 2011 when the trees were 13 month-old (trees were planted in November 2010).

  • Fundão (Portugal) in May 2013 when the trees were 27 month-old (trees were planted in February 2011)

3 Mortality models

In this section, we want to determine whether genomic offset (GO) or climate transfer distances (CTD) are associated with the proportion of dead trees in the populations, independently in the five common gardens. We build a model that assumes that the initial sapling height acts as a confounder. Indeed, trees that were higher at the time of planting have a higher probability of survival. This is particularly true in Madrid and Cáceres (Spain) where the seedlings experienced an extreme drought event the same year they were planted, resulting in very high mortality rates. Here is the model:

\[\begin{align*} a_{p} &\sim \text{Binomial} (N_{p},p_{p}) \\ \text{logit}(p_{p}) &= \beta_{0} + \beta_{H}H_{p} + \beta_{X}X_{p}\\ \end{align*}\]

with \(a_{p}\) the count of individuals that died in the population \(p\), \(N_{p}\) the total number of individuals in the population \(p\) (=number of individuals that were initially planted in the common garden), \(p_p\) is the estimated probability of mortality in the population \(p\), \(X_{p}\) is the genomic offset or climatic transfer distance for the population \(p\) and \(H_{p}\) is the BLUPs for height of the population \(p\) (population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. 2022). \(H_{p}\) is used as a proxy of the initial tree height.

We used the following weakly informative priors:

\[\begin{align*} \begin{bmatrix} \beta_{0,c} \\ \beta_{H} \\ \beta_{X} \end{bmatrix} &\sim \mathcal{N}(0,5) \end{align*}\]

3.1 The variables

3.1.1 Initial tree height (confounder)

We load the population-specific intercepts from the model 1 of Archambeau et al. (2022).

Code
mod_arch2022 <- readRDS(file=here("data/Archambeauetal2022_MOD1.rds"))
pop_heights <- mod_arch2022$fit %>% 
  broom.mixed::tidyMCMC(estimate.method = "mean",conf.int = T) %>% # we take the mean of the prov random intercepts
  filter(str_detect(term, "^(r_prov\\[)")) %>% 
  dplyr::rename(height=estimate,pop=term) %>% 
  mutate(pop=str_sub(pop,8,-12))

pop_heights[1:5,] %>% kable_mydf
pop height std.error conf.low conf.high
ALT 0.09 0.04 0.01 0.17
ARM 0.12 0.04 0.04 0.20
ARN -0.09 0.03 -0.16 -0.03
BAY -0.14 0.03 -0.20 -0.07
BON -0.04 0.04 -0.12 0.04

We export the height intercepts from Archambeau et al. (2022) to the DRYAD repository.

Code
pop_heights %>% write_csv(here("data/DryadRepo/HeightIntercepts_Archambeauetal2023.csv"))

3.1.2 Survival data

The response variable is counts of dead trees. To calculate these counts, we load the survival data in the common gardens, in which 0 corresponds to dead trees and 1 to survivors.

Code
surv_measurements <- c("AST_survmar14","BDX_surv18","CAC_survdec11","MAD_survdec11","POR_survmay13")
survival_data <- pheno_data %>% 
  dplyr::select(site,block,pop,clon,tree,any_of(surv_measurements)) %>% 
  pivot_longer(cols=any_of(surv_measurements), names_to = NULL, values_to = "survival") %>% 
  drop_na(survival) 

survival_data[1:5,] %>% kable_mydf
site block pop clon tree survival
asturias 2 CAR CAR12 CAR12_2 1
asturias 4 STJ STJ14 STJ14_4 1
asturias 5 CUE CUE6 CUE6_5 1
asturias 6 MAD MAD1 MAD1_6 1
asturias 4 ORI ORI12 ORI12_4 1

3.2 Run the models

Stan code of the mortality model:

Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_BinomialMortalityModel.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  int nb_tot[N];    // Total number of trees in the population
  int nb_dead[N];   // Number of dead trees in the population
  vector[N] H;      // Mean tree height of the population
  vector[N] X;      // Genomic offset or climatic transfer distance of the population
}
parameters {
  real beta_0;
  real beta_H;
  real beta_X;
}
model {
  nb_dead ~ binomial_logit(nb_tot,beta_0 + beta_H * H + beta_X * X);

  beta_0 ~ normal(0,5);//std_normal();
  beta_H ~ normal(0,5);//std_normal();
  beta_X ~ normal(0,5);//std_normal();
}
generated quantities{
  vector[N] log_lik;
  // log likelihood for loo
  for (n in 1:N) {
    log_lik[n] = binomial_logit_lpmf( nb_dead[n] | nb_tot[n] , beta_0 + beta_H * H[n] + beta_X * X[n]);
  }
} 

We run the models for each common garden and each genomic offset predictions/climatic transfer distances, so a total of 5 * 15 = 75 models runs.

Code
coefftab <- lapply(unique(survival_data$site),function(site_i){
  
  
  lapply(unique(df$method), function(method_i){
  
# Subset the data for the site i and method i
    sub_data <- survival_data %>% 
      filter(site == site_i) %>% 
      group_by(pop) %>% 
      dplyr::summarise(nb_dead=n()-sum(survival),nb_tot=n()) # transform survival data into mortality data
    
    sub_data <- df %>% 
      filter(method == method_i & cg == site_i) %>% 
      inner_join(sub_data, by="pop") %>% 
      inner_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))), by="pop") %>% 
      arrange(pop)
      
# Data in a list for Stan 
    stanlist <- list(N = nrow(sub_data),
                     nb_dead = sub_data$nb_dead,
                     nb_tot=sub_data$nb_tot,
                     H=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX))
    
# Running the model
    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE) 
    
    
    #loo(mod) %>% saveRDS(file=here(paste0("outputs/ValidationCommonGarden/MortalityModels/loos/loo_",site_i,"_",method_i,".rds")))
  
  
# Save coefficients
    broom.mixed::tidyMCMC(mod,
                  droppars = NULL, 
                  robust = FALSE, # give mean and standard deviation
                  ess = F, 
                  rhat = F, 
                  conf.int = T,
                  conf.level = 0.95) %>% 
    filter(str_detect(term, c('beta'))) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method=method_i,
           .before=1)
    
  }) %>% bind_rows()
  
}) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/MortalityModels/coefftab.rds"))

3.3 Model coefficients

We load the model coefficients:

Code
coefftab <- readRDS(file=here("outputs/ValidationCommonGarden/MortalityModels/coefftab.rds"))
coefftab[1:5,] %>% kable_mydf()
cg method term mean std_deviation conf_low conf_high
asturias CTD_bio1 beta_0 -2.99 0.08 -3.14 -2.84
asturias CTD_bio1 beta_H 0.15 0.07 0.01 0.28
asturias CTD_bio1 beta_X -0.02 0.07 -0.17 0.12
asturias CTD_bio12 beta_0 -2.98 0.07 -3.12 -2.84
asturias CTD_bio12 beta_H 0.14 0.07 -0.01 0.28

Below, we plot the 95% credible intervals of:

  • the \(\beta_X\) coefficients, which stand for the association between the counts of dead trees and the genomic offset predictions (i.e. capturing the potential maladaptation of the populations when transplanted in the new environment of the common gardens) or the climatic transfer distances.

  • the \(\beta_H\) coefficients, which stand for the association between the counts of dead trees and the BLUPs for height in the five common gardens, used as a proxy of the initial seedling height (a confounder in the model).

Code
coeff_match <- list(beta_H="$\\beta_{H}$ estimates (effect of initial tree height)",
                    beta_X="\\beta_{X}$ estimates") #$\\beta_{X}$ estimates
  
p <- lapply(c("beta_X","beta_H"), function(coeff){
  
p <- coefftab %>% 
  filter(term==coeff) %>% 
  left_join(distinct(df[,c("method_type","method_input","cg","method")]), by=c("method","cg")) %>% 
  mutate(cg=case_when(cg=="asturias"~"Asturias, Spain (37 months)",
                      cg=="bordeaux"~"Bordeaux, France (85 months)",
                      cg=="caceres"~"Cáceres, Spain (8 months)",
                      cg=="madrid"~"Madrid, Spain (13 months)",
                      cg=="portugal"~"Fundão, Portugal (27 months)"),
         method_input=case_when(method_input=="all" ~ "All SNPs",
                                method_input=="bio1" ~ "Mean annual temperature (°C)",
                                method_input=="bio12" ~ "Annual precipitation (mm)",
                                method_input=="bio15" ~ "Precipitation seasonality (index)",
                                method_input=="bio3" ~ "Isothermality (index)",
                                method_input=="bio4" ~ "Temperature seasonality (°C)",
                                method_input=="cand" ~ "Candidate SNPs",
                                method_input=="control" ~ "Control SNPs",
                                method_input=="SHM" ~ "Summer heat moisture index (°C/mm)")) %>% 
  mutate(method_input=factor(method_input, levels=c("All SNPs",
                                                    "Candidate SNPs",
                                                    "Control SNPs",
                                                    "Mean annual temperature (°C)",
                                                    "Annual precipitation (mm)",
                                                    "Precipitation seasonality (index)",
                                                    "Isothermality (index)",
                                                    "Temperature seasonality (°C)",
                                                    "Summer heat moisture index (°C/mm)")),
         method_type=factor(method_type, levels=c("GDM","GF","LFMM","RDA","CTD")))



# Plots with CTD and SNP sets
# ===========================
p_allvar <- p %>% ggplot(aes(x = method_type,
                             y = mean,
                             ymin = conf_low, 
                             ymax = conf_high,
                             color=method_input,
                             shape=method_input)) +
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = 0.6),
                     point_size=3.5, size=3) + # 
  facet_wrap(~cg, ncol=2) + #, scales="free_x", space = "free" 
  ylab(TeX(coeff_match[[coeff]])) + xlab("") +
  scale_color_manual(values=colors_coeff)+
  #paletteer::scale_color_paletteer_d("ggthemes::calc") +
  scale_shape_manual(values = c(rep(17,3),rep(16,6))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=12),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))

# save in pdf and png
ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.pdf")),
       device="pdf",
       height=9,
       width=12)

ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.png")),
       height=9,
       width=12)


# Plots with SNPs only
# ====================
p_SNPs <- p %>% 
  filter(!method_type == "CTD") %>% 
  ggplot(aes(x = method_type,
             y = mean,
             ymin = conf_low, 
             ymax = conf_high,
             color=method_input,
             shape=method_input)) +
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = 0.45),
                     point_size=3.5, size=3) + # 
  facet_wrap(~cg, ncol=2) + #, scales="free_x", space = "free" 
  ylab(TeX(coeff_match[[coeff]])) + xlab("") +
  scale_color_manual(values=colors_coeff)+
  #paletteer::scale_color_paletteer_d("ggthemes::calc") +
  scale_shape_manual(values = c(rep(17,3),rep(16,6))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))
  

# save in pdf and png
ggsave(p_SNPs,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsets.pdf")),
       device="pdf",
       height=7,
       width=8)

ggsave(p_SNPs,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsets.png")),
       height=7,
       width=8)


return(p_allvar)

})

p
[[1]]


[[2]]

Interpretation

As expected, the association between mortality rates and the initial tree height was particularly strong in Madrid and Cáceres (Spain) where the seedlings experienced an extreme drought event the same year they were planted, resulting in very high mortality rates. More surprisingly, this association was also osberved in Fundão, Portugal. However, the initial tree height was not associated with mortality rates in Bordeaux (France) and Asturias (Spain), which benefit from the favorable climates of the Atlantic region and in which the mortality rates were very low.

In Fundão, Cáceres and Madrid, for most genomic offset predictions, populations experiencing higher mortality rates also showed higher genomic offset. The most consistent associations across the three common gardens were obtained with:

  • Gradient Forest (GF) with both candidate and control SNPs.

  • Redundancy analysis (RDA) with both candidate and control SNPs.

  • Latent factor mixed model (LFMM) with all SNPs or control SNPs.

Interestingly, climatic transfer distances generally explain mortality probability less well than genomic offset predictions, and none of them show a consistent association with the number of dead trees across the three common gardens. Note that it is almost the case for the climatic transfer distance calculated based on the temperature seasonality, which shows a positive association with the counts of dead trees in Madrid and Fundão, and almost in Cáceres.

In Asturias, no association was detected between genomic offset predictions and mortality rates, which may be due to the favorable climatic conditions of this common garden and the associated low mortality rates (only 206 dead trees). Therefore, in this common garden, mortality events are likely to be mostly random events due to other factors than climate. However, we can note that an association was detected between mortality rates and isothermality: populations from areas with high isothermality tend to die more in the common garden of Asturias, Spain.

In Bordeaux, there were even less dead trees than in Asturias (119 dead trees). However, the two nonparametric approach used to predict the genomic offset (GDM and GF) detected a positive association between mortality rates and the genomic offset predictions obtained with both the candidate and the control SNPs. A positive association was also detected with the genomic offset predictions obtained with candidate SNPs with the LFMM approach. On the other hand, none of the genomic offset predictions from the RDA approach were associated with the mortality rates. Interestingly, the only association with a climatic transfer distance was with the temperature seasonality (bio4), which was the most important variable to explain the turnover in allele frequency in the GF and GDM approaches.

3.4 Experimental design

We export in a table with the number and proportion of dead trees for each population in each common garden.

Code
ExpDesignTab <- lapply(unique(survival_data$site),function(site_i){
  
survival_data %>% 
    filter(site==site_i) %>% 
    dplyr::select(pop,survival) %>% 
    drop_na() %>% 
    group_by(pop) %>% 
    dplyr::summarise(nb_dead=n()-sum(survival),nb_tot=n()) %>% 
    mutate(prop_dead=nb_dead*100/nb_tot)
  
}) %>% 
  setNames(unique(survival_data$site)) %>% 
  bind_rows(.id="cg") %>% 
  pivot_wider(names_from=cg,values_from = c(nb_dead, nb_tot,prop_dead),names_sep="_") %>% 
  dplyr::select(pop,contains("asturias"),contains("bordeaux"),contains("caceres"),contains("madrid"),contains("portugal"))

# To generate a latex table
# =========================
# print(xtable(ExpDesignTab, type = "latex",digits=2), 
#       file = paste0(here("tables/ExpDesignMortalityModelsPerPopCG.tex")), 
#       include.rownames=FALSE)

# Export the table for the Supplemetary Information
# =================================================
ExpDesignTab %>% 
  dplyr::select(-contains("prop_dead")) %>% 
  saveRDS(here("tables/ExpDesignMortalityModelsPerPopCG.rds"))


# Show table
# ==========
ExpDesignTab %>% kable_mydf()
pop nb_dead_asturias nb_tot_asturias prop_dead_asturias nb_dead_bordeaux nb_tot_bordeaux prop_dead_bordeaux nb_dead_caceres nb_tot_caceres prop_dead_caceres nb_dead_madrid nb_tot_madrid prop_dead_madrid nb_dead_portugal nb_tot_portugal prop_dead_portugal
ALT 3 72 4.17 1 53 1.89 69 72 95.83 52 72 72.22 19 72 26.39
ARM 1 64 1.56 1 64 1.56 56 64 87.50 49 64 76.56 23 64 35.94
ARN 5 136 3.68 3 109 2.75 129 136 94.85 107 136 78.68 46 136 33.82
BAY 8 144 5.56 3 142 2.11 132 144 91.67 113 144 78.47 56 144 38.89
BON 2 71 2.82 1 48 2.08 62 72 86.11 48 72 66.67 14 72 19.44
CAD 5 80 6.25 3 70 4.29 74 80 92.50 62 80 77.50 22 80 27.50
CAR 0 48 0.00 0 46 0.00 44 48 91.67 36 48 75.00 21 48 43.75
CAS 4 80 5.00 6 78 7.69 78 80 97.50 66 80 82.50 31 80 38.75
CEN 2 72 2.78 1 38 2.63 66 72 91.67 50 72 69.44 18 72 25.00
COC 6 144 4.17 3 109 2.75 137 144 95.14 118 144 81.94 53 143 37.06
COM 0 32 0.00 3 32 9.38 28 32 87.50 21 32 65.62 10 32 31.25
CUE 9 224 4.02 9 194 4.64 211 224 94.20 186 224 83.04 101 224 45.09
HOU 11 208 5.29 6 160 3.75 186 208 89.42 145 208 69.71 66 208 31.73
LAM 3 72 4.17 4 60 6.67 68 72 94.44 57 72 79.17 27 72 37.50
LEI 14 184 7.61 10 143 6.99 162 184 88.04 147 184 79.89 64 184 34.78
MAD 0 8 0.00 0 8 0.00 8 8 100.00 5 8 62.50 2 8 25.00
MIM 7 144 4.86 4 126 3.17 135 144 93.75 116 144 80.56 67 144 46.53
OLB 3 176 1.70 4 116 3.45 149 176 84.66 113 176 64.20 27 176 15.34
OLO 22 191 11.52 3 132 2.27 173 192 90.10 145 192 75.52 73 192 38.02
ORI 4 208 1.92 2 183 1.09 196 208 94.23 162 208 77.88 67 208 32.21
PET 11 192 5.73 3 147 2.04 171 192 89.06 136 192 70.83 75 192 39.06
PIA 8 128 6.25 1 110 0.91 111 128 86.72 91 128 71.09 43 128 33.59
PIE 3 72 4.17 4 55 7.27 67 72 93.06 50 72 69.44 17 72 23.61
PLE 8 160 5.00 9 138 6.52 155 160 96.88 131 160 81.88 71 160 44.38
PUE 1 64 1.56 5 56 8.93 58 64 90.62 50 64 78.12 30 64 46.88
QUA 5 136 3.68 3 117 2.56 123 136 90.44 83 136 61.03 36 136 26.47
SAC 1 72 1.39 2 51 3.92 65 72 90.28 56 72 77.78 35 72 48.61
SAL 7 112 6.25 0 69 0.00 104 112 92.86 85 112 75.89 55 112 49.11
SEG 12 168 7.14 7 168 4.17 154 168 91.67 135 168 80.36 77 168 45.83
SIE 2 64 3.12 1 47 2.13 59 64 92.19 40 64 62.50 29 64 45.31
STJ 15 224 6.70 3 180 1.67 190 224 84.82 154 224 68.75 70 224 31.25
TAM 8 120 6.67 4 71 5.63 114 120 95.00 107 120 89.17 60 120 50.00
VAL 5 96 5.21 5 64 7.81 87 96 90.62 69 96 71.88 29 96 30.21
VER 11 216 5.09 5 160 3.12 197 216 91.20 153 216 70.83 83 216 38.43
Code
# Information used in the manuscript
ExpDesignTab[,-1] %>% 
  dplyr::summarise_all(mean) %>% 
  kable_mydf
nb_dead_asturias nb_tot_asturias prop_dead_asturias nb_dead_bordeaux nb_tot_bordeaux prop_dead_bordeaux nb_dead_caceres nb_tot_caceres prop_dead_caceres nb_dead_madrid nb_tot_madrid prop_dead_madrid nb_dead_portugal nb_tot_portugal prop_dead_portugal
6.06 123 4.27 3.5 98.35 3.7 112.29 123.06 91.65 92.29 123.06 74.31 44.62 123.03 35.79

4 Height models

In this section, we estimate the association between genomic offset (GO) predictions or climate transfer distances (CTD) and tree height, independently in the five common gardens. We compare four different mathematical models.

We first subset height measurements from the dataset.

Code
height_measurements <- c("AST_htmar14","BDX_htnov18","CAC_htdec11","MAD_htdec11","POR_htmay13")

height_data <- pheno_data %>% 
  dplyr::rename(cg = site) %>% 
  dplyr::select(cg,block,pop,clon,tree,any_of(height_measurements)) %>% 
  pivot_longer(cols=any_of(height_measurements), names_to=NULL,values_to="height",values_drop_na = TRUE)

height_data[1:5,] %>% kable_mydf
cg block pop clon tree height
asturias 2 CAR CAR12 CAR12_2 840
asturias 4 STJ STJ14 STJ14_4 1200
asturias 5 CUE CUE6 CUE6_5 1160
asturias 6 MAD MAD1 MAD1_6 1340
asturias 4 ORI ORI12 ORI12_4 1320

We write a function to plot the model coefficients :

Code
coeff_y_axis <- list(beta_X1="Regression coefficients $\\beta_{X_1}$",
                     beta_X2="Regression coefficients $\\beta_{X_2}$",
                     beta_H="Regression coefficients $\\beta_H$")

generate_interval_plots <- function(model_i){

coefftab <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/coefftab_",model_i,".rds")))

  
params <-  coefftab$term %>% unique() %>% str_subset("^beta")

p <- lapply(params, function(coeff){


p <- coefftab %>% 
  filter(term==coeff) %>% 
  left_join(distinct(df[,c("method_type","method_input","cg","method")]), by=c("method","cg")) %>% 
  mutate(cg=case_when(cg=="asturias"~"Asturias, Spain (37 months)",
                      cg=="bordeaux"~"Bordeaux, France (85 months)",
                      cg=="caceres"~"Cáceres, Spain (8 months)",
                      cg=="madrid"~"Madrid, Spain (13 months)",
                      cg=="portugal"~"Fundão, Portugal (27 months)"),
         method_input=case_when(method_input=="all" ~ "All SNPs",
                                method_input=="bio1" ~ "Mean annual temperature (°C)",
                                method_input=="bio12" ~ "Annual precipitation (mm)",
                                method_input=="bio15" ~ "Precipitation seasonality (index)",
                                method_input=="bio3" ~ "Isothermality (index)",
                                method_input=="bio4" ~ "Temperature seasonality (°C)",
                                method_input=="cand" ~ "Candidate SNPs",
                                method_input=="control" ~ "Control SNPs",
                                method_input=="SHM" ~ "Summer heat moisture index (°C/mm)")) %>% 
  mutate(method_input=factor(method_input, levels=c("All SNPs",
                                                    "Candidate SNPs",
                                                    "Control SNPs",
                                                    "Mean annual temperature (°C)",
                                                    "Annual precipitation (mm)",
                                                    "Precipitation seasonality (index)",
                                                    "Isothermality (index)",
                                                    "Temperature seasonality (°C)",
                                                    "Summer heat moisture index (°C/mm)")),
         method_type=factor(method_type, levels=c("GDM","GF","LFMM","RDA","CTD")))
  

# Plots with CTD and SNP sets
# ===========================
p_allvar <- p %>% 
   ggplot(aes(x = method_type,
              y = mean,
              ymin = conf_low, 
              ymax = conf_high,
              color=method_input,
              shape=method_input)) +
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) + # 
  facet_wrap(~cg, ncol=2) + #, scales="free_x", space = "free" 
  ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
  scale_color_manual(values=colors_coeff) +
  scale_shape_manual(values = c(rep(16,6),rep(17,3))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=12),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))


# save in pdf
ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD.pdf")),
       device="pdf",
       height=9,
       width=12)

# save in png
# ggsave(p_allvar,
#        file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD.png")),
#        height=9,
#        width=12)


# Plots only SNP sets
# ===================
p_SNPs <- p %>% 
  filter(!method_type == "CTD") %>% 
  ggplot(aes(x = method_type,
             y = mean,
             ymin = conf_low, 
             ymax = conf_high,
             color=method_input,
             shape=method_input)) +
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) + # 
  facet_wrap(~cg, ncol=2) + #, scales="free_x", space = "free" 
  ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
  scale_color_manual(values=colors_coeff) +
  scale_shape_manual(values = c(rep(16,6),rep(17,3))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))


# save in pdf
ggsave(p_SNPs,
       file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsets.pdf")),
       device="pdf",
       height=7,
       width=8)

# save in png
# ggsave(p_SNPs,
#        file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsets.png")),
#        height=7,
#        width=8)

return(p_allvar)

})

return(p)
}

We write a function to plot the quadratic relationships between tree height and the GO predictions / CTD.

Code
generate_poly_plots <- function(model_i){

coefftab <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/coefftab_",model_i,".rds")))

x <- seq(-2,2,0.01)
poly_function <- function(a,b,x) {a * x^2 + b * x}

ggtab <- coefftab %>% 
  filter(term %in% c("beta_X1","beta_X2")) %>% 
  left_join(distinct(df[,c("method_type","method_input","cg","method")]), by=c("method","cg")) %>% 
  mutate(cg_legend =case_when(cg=="asturias"~"Asturias (37 months)",
                      cg=="bordeaux"~"Bordeaux (85 months)",
                      cg=="caceres"~"Cáceres (8 months)",
                      cg=="madrid"~"Madrid (13 months)",
                      cg=="portugal"~"Fundão (27 months)"),
         method_input=case_when(method_input=="all" ~ "All SNPs",
                                method_input=="bio1" ~ "Mean annual temperature (°C)",
                                method_input=="bio12" ~ "Annual precipitation (mm)",
                                method_input=="bio15" ~ "Precipitation seasonality (index)",
                                method_input=="bio3" ~ "Isothermality (index)",
                                method_input=="bio4" ~ "Temperature seasonality (°C)",
                                method_input=="cand" ~ "Candidate SNPs",
                                method_input=="control" ~ "Control SNPs",
                                method_input=="SHM" ~ "Summer heat moisture index (°C/mm)")) %>% 
  pivot_wider(names_from="term", values_from = c("mean","std_deviation","conf_low","conf_high"))



ggtab <- lapply(unique(ggtab$cg), function(cg_i){
lapply(unique(ggtab$method), function(method_i){
  
  subggtab <- ggtab %>% filter(method == method_i & cg == cg_i)
  
  poly_function(a=subggtab$mean_beta_X2, b=subggtab$mean_beta_X1,x=x)
  }) %>% 
    setNames(unique(ggtab$method)) %>% 
    bind_rows(.id="method") %>% 
    mutate(x_axis = x)
}) %>% 
  setNames(unique(ggtab$cg)) %>% 
  bind_rows(.id="cg") %>% 
  pivot_longer(cols=any_of(unique(ggtab$method)),names_to="method",values_to="predictions") %>% 
  left_join(ggtab[,c("method","method_type", "method_input")] %>% distinct(), by="method")
  

# Below, you can use unique(ggtab$cg) or unique(ggtab$cg_legend)
ggplots <- lapply(unique(ggtab$cg), function(cg_i){
  
p_go <- ggtab %>%
  filter(cg %in% cg_i & method_type %in% unique(ggtab$method_type)[!unique(ggtab$method_type) == "CTD"]) %>% 
  mutate(method_input = factor(method_input, levels=c("Candidate SNPs", "Control SNPs", "All SNPs"))) %>% 
  ggplot(aes(x=x_axis,y=predictions,col=method_type,linetype=method_input)) +
  geom_hline(yintercept=0,color="black") +
  geom_line(linewidth=1) +
  ylab('Standardized tree height') + 
  #ylab(paste0('Standardized tree height - ',cg_i)) + 
  xlab(TeX("Standardized predictor")) +
  ylim(c(min(ggtab$predictions),max(ggtab$predictions))) +
  scale_color_manual(values=c("#FCA315FF", "#00FF00FF","#172869FF", "#FF00FFFF"), name = "GEA method") +
  scale_linetype_manual(values=c("solid","dashed","dotted"), name = "SNP set") +
  theme_bw() + 
  theme(axis.text.x = element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title = element_text(size=16),
        legend.box="horizontal",
        legend.key.width = unit(1,"cm"),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        legend.position = c(0.73,0.89)) + 
  guides(color=guide_legend(ncol=1))

#if(!cg_i==unique(ggtab$cg)[1])  p_go <- p_go + theme(legend.position = "none")
#if(!cg_i==unique(ggtab$cg)[length(unique(ggtab$cg))])  p_go <- p_go + xlab("")

p_ctd <- ggtab %>% 
  filter(cg %in% cg_i & method_type %in% "CTD") %>%
  ggplot(aes(x=x_axis,y=predictions,col=method_input)) +
  geom_hline(yintercept=0,color="black") +
  #geom_segment(aes(x = min(x_axis), y = 0, xend = max(x_axis), yend = 0), color="black") +
  geom_line(linewidth=1.2) +
  ylab('') + 
  xlab(TeX("Standardized predictor")) +
  ylim(c(min(ggtab$predictions),max(ggtab$predictions))) +
  scale_color_manual(values=colors_coeff[4:length((colors_coeff))], 
                     name = "Climatic transfer distance") +
    theme_bw() + 
    theme(axis.text.x = element_text(size=13),
          axis.text.y = element_text(size=13),
          axis.title = element_text(size=16),
          legend.box="horizontal",
          legend.background = element_blank(),
          legend.box.background = element_blank(),
          legend.key = element_blank(),
          legend.position = c(0.55,0.91)) + 
  guides(color=guide_legend(ncol=2))

#if(!cg_i==unique(ggtab$cg)[1])  p_ctd <- p_ctd + theme(legend.position = "none")
#if(!cg_i==unique(ggtab$cg)[length(unique(ggtab$cg))])  p_ctd <- p_ctd + xlab("")

pp <- plot_grid(p_go,p_ctd)

ggsave(filename = here(paste0("figs/ValidationCommonGarden/HeightModels/ScatterPlotsPredictedRelationships_",model_i,"_",cg_i,".pdf")), device="pdf", width=12,height=7)

pp

})

return(ggplots)
}

4.1 Model 1 (basic one)

4.1.1 Model equation and code

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= B_b + \beta_{X1}X_p + \beta_{X2}X^2_p\\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M1.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;                                    // Number of individuals
  vector[N] Y;                              // Response variable (individual tree height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real<lower = 0>  sigma;                   // Residual variance of the model
}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model

  mu = alpha_bloc[bloc] + beta_X1 * X + beta_X2 * square(X);
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","R_squared","sigma","alpha_bloc")
Code
coefftab <- lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method), function(method_i){
  
    df_sub <- df %>% filter(method == method_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop"))
      
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     nb_bloc = length(unique(sub_data$block)),
                     bloc = as.numeric(as.factor(sub_data$block)))

    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE,
                    pars=params_to_estimate) 
  
# Model coefficients
    broom.mixed::tidyMCMC(mod,
                          pars=params_to_estimate,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method=method_i,
           .before=1)
    
    
  }) %>% bind_rows()
  
}) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M1.rds"))

4.1.2 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\) and \(\beta_{X_2}\) coefficients.

Code
generate_interval_plots(model_i = "M1")
[[1]]


[[2]]

4.1.3 Predicted quadratic relationships

We plot the predicted quadratic relationships between tree height and GO predictions / CTD.

Code
generate_poly_plots(model_i="M1")
[[1]]


[[2]]


[[3]]


[[4]]


[[5]]

4.1.4 Interpretation

Interpretation

Higher genomic offset predictions are consistently associated with lower tree height only in Asturias. Most genomic offset predictions are also associated with lower tree height in Bordeaux.

In Fundão, Cáceres and Madrid, the association between genomic offset predictions and tree height are often in the opposite direction as expected.

Interestingly, the climatic transfer distances based in the mean annual temperature and the annual precipitation are negatively associated with tree height in all common gardens (except Cáceres for the CTD based on annual precipitation). These CTD may therefore be better predictors of tree height in common gardens than the genomic offset predictions.

Finally, I think that using tree height as a proxy for fitness to evaluate genomic offset predictions in common gardens may not be appropriate in maritime pine because this trait has a very low genetic-by-environment interaction (see previous papers). Therefore, the association between genomic offset predictions and tree height in Asturias and Bordeaux may not be due to a higher fitness of the trees that are best adapted to the climate in these common gardens, but more probably to height differences that are common across all environments (i.e. populations from Atlantic climates are generally taller). This is confirmed by the models below in which the population-specific BLUPs for height across all common gardens was included as a confounder. When included, the association between genomic offset predictions and tree height in Asturias almost disappears and either disappears or goes in the opposite direction as expected in Bordeaux. Except the genomic offset predictions based on the RDA remain negatively associated with tree height in Bordeaux and Asturias.

Therefore, I am not sure we can go very far in this validation step. We may retain that genomic offset predictions seem to show the most consistent association with tree height when generated with the RDA, and that this association is robust (though smaller) even when the fixed genetic differences across populations are included as confounder.

4.2 Model 2

4.2.1 Model equation and code

Model 2 = Model 1 + Initial height as a confounder.

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M2.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model

}
transformed parameters {
  vector[N] mu;    // Linear predictor
  real R_squared;  // R^2 to evaluate the goodness of fit of the model

  mu = alpha_bloc[bloc] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  beta_H ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","R_squared","sigma","alpha_bloc")
Code
coefftab <- lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method), function(method_i){
  
    df_sub <- df %>% filter(method == method_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
      
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     bloc = as.numeric(as.factor(sub_data$block)))

    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE,
                    pars = params_to_estimate) 
  
# Save coefficients
    broom.mixed::tidyMCMC(mod,
                          pars=params_to_estimate,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method=method_i,
           .before=1)
    
    
  }) %>% bind_rows()
  
}) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M2.rds"))

4.2.2 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients.

Code
generate_interval_plots(model_i = "M2")
[[1]]


[[2]]


[[3]]

4.2.3 Predicted quadratic relationships

Code
generate_poly_plots(model_i = "M2")
[[1]]


[[2]]


[[3]]


[[4]]


[[5]]

4.3 Model 3

4.3.1 Model equation and code

Model 3 = Model 2 + Initial height as a confounder + Clone fixed intercepts

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= \beta_0 + B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \beta_0 &\sim \mathcal{N}(\text{mean}(Y),2) \\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ C_c \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
  • \(\beta_0\) is the model global intercept.
  • \(C_c\) are the clone intercepts.
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M3.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
  int<lower=0> nb_clon;                     // Number of clones
  int<lower=0, upper=nb_clon> clon[N];      // Clones
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  vector[nb_clon] alpha_clon;               // Intercepts of the clones
  real beta_0;                              // Global intercept
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Linear coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model
}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model

  mu = beta_0 + alpha_bloc[bloc] + alpha_clon[clon] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  alpha_clon ~ std_normal();
  beta_0 ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
  beta_H ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","R_squared","sigma","alpha_bloc","alpha_clon")
Code
coefftab <- lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method), function(method_i){
  
    df_sub <- df %>% filter(method == method_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")

# Datalist for Stan      
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     nb_clon = length(unique(sub_data$clon)),
                     bloc = as.numeric(as.factor(sub_data$block)),
                     clon = as.numeric(as.factor(sub_data$clon)))
    
# Running the model
    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE,
                    pars = params_to_estimate) 
  
# Model coefficients
    broom.mixed::tidyMCMC(mod,
                          pars=params_to_estimate,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method=method_i,
           .before=1)
    
    
  }) %>% bind_rows()
  
}) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M3.rds"))

4.3.2 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients.

Code
generate_interval_plots(model_i = "M3")
[[1]]


[[2]]


[[3]]

4.3.3 Predicted quadratic relationships

Code
generate_poly_plots(model_i = "M3")
[[1]]


[[2]]


[[3]]


[[4]]


[[5]]

4.4 Model 4

4.4.1 Model equation and code

Model 4 = Model 2 + Initial height as a confounder + Clone varying intercepts

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^2) \\ \mu_{pb} &= \beta_0 + B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \beta_0 &\sim \mathcal{N}(\text{mean}(Y),2) \\ C_c &\sim \mathcal{N}(0,\sigma_C^2) \\ \begin{bmatrix} \sigma \\ \sigma_C \end{bmatrix} &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
  • \(\beta_0\) is the model global intercept.
  • \(C_c\) are the clone varying intercepts which follow a normal distribution of variance \(\sigma_C^2\).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M4.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
  int<lower=0> nb_clon;                     // Number of clones
  int<lower=0, upper=nb_clon> clon[N];      // Clones
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_0;                              // Global intercept
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model
  real<lower = 0> sigma_clon;               // Variance of the clone intercepts
  vector[nb_clon] z_clon;                     // Vector for non-centered parameterization

}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model
  vector[nb_clon] alpha_clon;               // Intercepts of the clones

  alpha_clon = z_clon*sigma_clon;
  mu = beta_0 + alpha_bloc[bloc] + alpha_clon[clon] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  beta_0 ~ normal(mean(Y),2);
  sigma ~ exponential(1);
  sigma_clon ~ exponential(1);
  alpha_bloc ~ std_normal();
  z_clon ~ std_normal();
  beta_H ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","R_squared","sigma","alpha_bloc","alpha_clon")
Code
coefftab <- lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method), function(method_i){
  
    df_sub <- df %>% filter(method == method_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")

# Datalist for Stan      
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     nb_clon = length(unique(sub_data$clon)),
                     bloc = as.numeric(as.factor(sub_data$block)),
                     clon = as.numeric(as.factor(sub_data$clon)))
    
# Running the model
    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE,
                    pars = params_to_estimate) 
  
# Model coefficients
    broom.mixed::tidyMCMC(mod,
                          pars=params_to_estimate,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method=method_i,
           .before=1)
    
    
  }) %>% bind_rows()
  
}) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M4.rds"))

4.4.2 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients.

Code
generate_interval_plots(model_i = "M4")
[[1]]


[[2]]


[[3]]

4.4.3 Predicted quadratic relationships

Code
generate_poly_plots(model_i = "M4")
[[1]]


[[2]]


[[3]]


[[4]]


[[5]]

4.5 Sample size for each clone

We look at the number of height measurements for each clone. In the tables below (one for each common garden), the first column correspond to the number of height measurements and the second column correspond to the number of clones with this number of measurements.

Code
lapply(unique(height_data$cg), function(site_i){
  
  
height_data %>% 
    filter(cg==site_i) %>% 
    dplyr::select(pop,clon) %>% 
    drop_na() %>% 
    group_by(clon) %>% 
    dplyr::summarise(nb_measurements=n()) %>% 
    count(nb_measurements)
  
  
}) %>% set_names(unique(height_data$cg))
$asturias
# A tibble: 5 × 2
  nb_measurements     n
            <int> <int>
1               4     1
2               5     5
3               6    30
4               7   129
5               8   358

$bordeaux
# A tibble: 5 × 2
  nb_measurements     n
            <int> <int>
1               4     3
2               5    13
3               6    36
4               7   108
5               8   271

$caceres
# A tibble: 4 × 2
  nb_measurements     n
            <int> <int>
1               1   183
2               2    58
3               3    11
4               4     2

$madrid
# A tibble: 6 × 2
  nb_measurements     n
            <int> <int>
1               1   148
2               2   149
3               3   101
4               4    51
5               5    15
6               6     3

$portugal
# A tibble: 8 × 2
  nb_measurements     n
            <int> <int>
1               1     7
2               2    20
3               3    49
4               4    91
5               5   135
6               6   127
7               7    66
8               8    26

5 DRYAD repository

We export mortality and height data in the DRYAD repository.

Code
pheno_data %>% 
  dplyr::rename(cg = site) %>% 
  dplyr::select(cg,block,pop,clon,tree,any_of(height_measurements),any_of(surv_measurements)) %>% 
  write_csv(here("data/DryadRepo/CommonGardenData_cleaned.csv"))

6 Session information

Code
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       Ubuntu 22.04.3 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2024-01-18
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version  date (UTC) lib source
 abind            1.4-5    2016-07-21 [1] CRAN (R 4.3.2)
 arrayhelpers     1.1-0    2020-02-04 [1] CRAN (R 4.3.2)
 backports        1.4.1    2021-12-13 [1] CRAN (R 4.3.2)
 bayesplot      * 1.10.0   2022-11-16 [1] CRAN (R 4.3.2)
 bit              4.0.5    2022-11-15 [1] CRAN (R 4.3.2)
 bit64            4.0.5    2020-08-30 [1] CRAN (R 4.3.2)
 broom            1.0.5    2023-06-09 [1] CRAN (R 4.3.2)
 broom.mixed      0.2.9.4  2022-04-17 [1] CRAN (R 4.3.2)
 cachem           1.0.8    2023-05-01 [1] CRAN (R 4.3.2)
 cellranger       1.1.0    2016-07-27 [1] CRAN (R 4.3.2)
 checkmate        2.3.1    2023-12-04 [1] CRAN (R 4.3.2)
 class            7.3-22   2023-05-03 [4] CRAN (R 4.3.1)
 classInt         0.4-10   2023-09-05 [1] CRAN (R 4.3.2)
 cli              3.6.2    2023-12-11 [1] CRAN (R 4.3.2)
 coda             0.19-4   2020-09-30 [1] CRAN (R 4.3.2)
 codetools        0.2-19   2023-02-01 [4] CRAN (R 4.2.2)
 colorspace       2.1-0    2023-01-23 [1] CRAN (R 4.3.2)
 cowplot        * 1.1.2    2023-12-15 [1] CRAN (R 4.3.2)
 crayon           1.5.2    2022-09-29 [1] CRAN (R 4.3.2)
 DBI              1.2.0    2023-12-21 [1] CRAN (R 4.3.2)
 devtools         2.4.5    2022-10-11 [1] CRAN (R 4.3.2)
 digest           0.6.33   2023-07-07 [1] CRAN (R 4.3.2)
 distributional   0.3.2    2023-03-22 [1] CRAN (R 4.3.2)
 dplyr          * 1.1.4    2023-11-17 [1] CRAN (R 4.3.2)
 e1071            1.7-14   2023-12-06 [1] CRAN (R 4.3.2)
 ellipsis         0.3.2    2021-04-29 [1] CRAN (R 4.3.2)
 evaluate         0.23     2023-11-01 [1] CRAN (R 4.3.2)
 fansi            1.0.6    2023-12-08 [1] CRAN (R 4.3.2)
 farver           2.1.1    2022-07-06 [1] CRAN (R 4.3.2)
 fastmap          1.1.1    2023-02-24 [1] CRAN (R 4.3.2)
 forcats        * 1.0.0    2023-01-29 [1] CRAN (R 4.3.2)
 fs               1.6.3    2023-07-20 [1] CRAN (R 4.3.2)
 furrr            0.3.1    2022-08-15 [1] CRAN (R 4.3.2)
 future           1.33.1   2023-12-22 [1] CRAN (R 4.3.2)
 generics         0.1.3    2022-07-05 [1] CRAN (R 4.3.2)
 ggdist           3.3.1    2023-11-27 [1] CRAN (R 4.3.2)
 ggplot2        * 3.4.4    2023-10-12 [1] CRAN (R 4.3.2)
 globals          0.16.2   2022-11-21 [1] CRAN (R 4.3.2)
 glue             1.6.2    2022-02-24 [1] CRAN (R 4.3.2)
 gridExtra        2.3      2017-09-09 [1] CRAN (R 4.3.2)
 gtable           0.3.4    2023-08-21 [1] CRAN (R 4.3.2)
 here           * 1.0.1    2020-12-13 [1] CRAN (R 4.3.2)
 highr            0.10     2022-12-22 [1] CRAN (R 4.3.2)
 hms              1.1.3    2023-03-21 [1] CRAN (R 4.3.2)
 htmltools        0.5.7    2023-11-03 [1] CRAN (R 4.3.2)
 htmlwidgets      1.6.4    2023-12-06 [1] CRAN (R 4.3.2)
 httpuv           1.6.13   2023-12-06 [1] CRAN (R 4.3.2)
 httr             1.4.7    2023-08-15 [1] CRAN (R 4.3.2)
 inline           0.3.19   2021-05-31 [1] CRAN (R 4.3.2)
 jsonlite         1.8.8    2023-12-04 [1] CRAN (R 4.3.2)
 kableExtra     * 1.3.4    2021-02-20 [1] CRAN (R 4.3.2)
 KernSmooth       2.23-22  2023-07-10 [4] CRAN (R 4.3.1)
 knitr          * 1.45     2023-10-30 [1] CRAN (R 4.3.2)
 labeling         0.4.3    2023-08-29 [1] CRAN (R 4.3.2)
 later            1.3.2    2023-12-06 [1] CRAN (R 4.3.2)
 latex2exp      * 0.9.6    2022-11-28 [1] CRAN (R 4.3.2)
 lattice          0.22-5   2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle        1.0.4    2023-11-07 [1] CRAN (R 4.3.2)
 listenv          0.9.0    2022-12-16 [1] CRAN (R 4.3.2)
 loo              2.6.0    2023-03-31 [1] CRAN (R 4.3.2)
 lubridate      * 1.9.3    2023-09-27 [1] CRAN (R 4.3.2)
 magrittr       * 2.0.3    2022-03-30 [1] CRAN (R 4.3.2)
 matrixStats      1.2.0    2023-12-11 [1] CRAN (R 4.3.2)
 memoise          2.0.1    2021-11-26 [1] CRAN (R 4.3.2)
 mime             0.12     2021-09-28 [1] CRAN (R 4.3.2)
 miniUI           0.1.1.1  2018-05-18 [1] CRAN (R 4.3.2)
 munsell          0.5.0    2018-06-12 [1] CRAN (R 4.3.2)
 nlme             3.1-163  2023-08-09 [4] CRAN (R 4.3.1)
 paletteer      * 1.5.0    2022-10-19 [1] CRAN (R 4.3.2)
 parallelly       1.36.0   2023-05-26 [1] CRAN (R 4.3.2)
 pillar           1.9.0    2023-03-22 [1] CRAN (R 4.3.2)
 pkgbuild         1.4.3    2023-12-10 [1] CRAN (R 4.3.2)
 pkgconfig        2.0.3    2019-09-22 [1] CRAN (R 4.3.2)
 pkgload          1.3.3    2023-09-22 [1] CRAN (R 4.3.2)
 plyr             1.8.9    2023-10-02 [1] CRAN (R 4.3.2)
 posterior        1.5.0    2023-10-31 [1] CRAN (R 4.3.2)
 profvis          0.3.8    2023-05-02 [1] CRAN (R 4.3.2)
 promises         1.2.1    2023-08-10 [1] CRAN (R 4.3.2)
 proxy            0.4-27   2022-06-09 [1] CRAN (R 4.3.2)
 purrr          * 1.0.2    2023-08-10 [1] CRAN (R 4.3.2)
 QuickJSR         1.0.9    2023-12-18 [1] CRAN (R 4.3.2)
 R6               2.5.1    2021-08-19 [1] CRAN (R 4.3.2)
 ragg             1.2.7    2023-12-11 [1] CRAN (R 4.3.2)
 raster         * 3.6-26   2023-10-14 [1] CRAN (R 4.3.2)
 RColorBrewer   * 1.1-3    2022-04-03 [1] CRAN (R 4.3.2)
 Rcpp             1.0.11   2023-07-06 [1] CRAN (R 4.3.2)
 RcppParallel     5.1.7    2023-02-27 [1] CRAN (R 4.3.2)
 readr          * 2.1.4    2023-02-10 [1] CRAN (R 4.3.2)
 readxl         * 1.4.3    2023-07-06 [1] CRAN (R 4.3.2)
 rematch2         2.1.2    2020-05-01 [1] CRAN (R 4.3.2)
 remotes          2.4.2.1  2023-07-18 [1] CRAN (R 4.3.2)
 reshape2       * 1.4.4    2020-04-09 [1] CRAN (R 4.3.2)
 rlang            1.1.2    2023-11-04 [1] CRAN (R 4.3.2)
 rmarkdown        2.25     2023-09-18 [1] CRAN (R 4.3.2)
 rnaturalearth  * 1.0.1    2023-12-15 [1] CRAN (R 4.3.2)
 rprojroot        2.0.4    2023-11-05 [1] CRAN (R 4.3.2)
 rstan          * 2.32.3   2023-10-15 [1] CRAN (R 4.3.2)
 rstudioapi       0.15.0   2023-07-07 [1] CRAN (R 4.3.2)
 rvest            1.0.3    2022-08-19 [1] CRAN (R 4.3.2)
 scales           1.3.0    2023-11-28 [1] CRAN (R 4.3.2)
 sessioninfo      1.2.2    2021-12-06 [1] CRAN (R 4.3.2)
 sf               1.0-15   2023-12-18 [1] CRAN (R 4.3.2)
 shiny            1.8.0    2023-11-17 [1] CRAN (R 4.3.2)
 sp             * 2.1-2    2023-11-26 [1] CRAN (R 4.3.2)
 StanHeaders    * 2.26.28  2023-09-07 [1] CRAN (R 4.3.2)
 stringi          1.8.3    2023-12-11 [1] CRAN (R 4.3.2)
 stringr        * 1.5.1    2023-11-14 [1] CRAN (R 4.3.2)
 svglite          2.1.3    2023-12-08 [1] CRAN (R 4.3.2)
 svUnit           1.0.6    2021-04-19 [1] CRAN (R 4.3.2)
 systemfonts      1.0.5    2023-10-09 [1] CRAN (R 4.3.2)
 tensorA          0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
 terra            1.7-65   2023-12-15 [1] CRAN (R 4.3.2)
 textshaping      0.3.7    2023-10-09 [1] CRAN (R 4.3.2)
 tibble         * 3.2.1    2023-03-20 [1] CRAN (R 4.3.2)
 tidybayes      * 3.0.6    2023-08-12 [1] CRAN (R 4.3.2)
 tidyr          * 1.3.0    2023-01-24 [1] CRAN (R 4.3.2)
 tidyselect       1.2.0    2022-10-10 [1] CRAN (R 4.3.2)
 tidyverse      * 2.0.0    2023-02-22 [1] CRAN (R 4.3.2)
 timechange       0.2.0    2023-01-11 [1] CRAN (R 4.3.2)
 tzdb             0.4.0    2023-05-12 [1] CRAN (R 4.3.2)
 units            0.8-5    2023-11-28 [1] CRAN (R 4.3.2)
 urlchecker       1.0.1    2021-11-30 [1] CRAN (R 4.3.2)
 usethis          2.2.2    2023-07-06 [1] CRAN (R 4.3.2)
 utf8             1.2.4    2023-10-22 [1] CRAN (R 4.3.2)
 vctrs            0.6.5    2023-12-01 [1] CRAN (R 4.3.2)
 viridisLite      0.4.2    2023-05-02 [1] CRAN (R 4.3.2)
 vroom            1.6.5    2023-12-05 [1] CRAN (R 4.3.2)
 webshot          0.5.5    2023-06-26 [1] CRAN (R 4.3.2)
 withr            2.5.2    2023-10-30 [1] CRAN (R 4.3.2)
 xfun             0.41     2023-11-01 [1] CRAN (R 4.3.2)
 xml2             1.3.6    2023-12-04 [1] CRAN (R 4.3.2)
 xtable         * 1.8-4    2019-04-21 [1] CRAN (R 4.3.2)
 yaml             2.3.8    2023-12-11 [1] CRAN (R 4.3.2)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Archambeau, Juliette, Marta Benito Garzón, Frédéric Barraquand, Marina de Miguel, Christophe Plomion, and Santiago C González-Martı́nez. 2022. “Combining Climatic and Genomic Data Improves Range-Wide Tree Height Growth Prediction in a Forest Tree.” The American Naturalist 200 (4): E141–59.
Fitzpatrick, Matthew C, Vikram E Chhatre, Raju Y Soolanayakanahally, and Stephen R Keller. 2021. “Experimental Support for Genomic Prediction of Climate Maladaptation Using the Machine Learning Approach Gradient Forests.” Molecular Ecology Resources.